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ABSTRACT 

We consider the early stages of cosmic hydrogen or helium reionization, when ionizing 
sources were still rare. We show that Poisson fluctuations in the galaxy distribution 
substantially affected the early bubble size distribution, although galaxy clustering was 
also an essential factor even at the earliest times. We find that even at high redshifts, 
a significant fraction of the ionized volume resided in bubbles containing multiple 
sources, regardless of the ionizing efficiency of sources or of the reionization redshift. In 
particular, for helium reionization by quasars, one-source bubbles last dominated (i.e., 
contained 90% of the ionized volume) at some redshift above z — 7.3, and hydrogen 
reionization by stars achieved this milestone at z > 23. For the early generations of 
atomic-cooling halos or molecular-hydrogen-cooling halos, one-source ionized regions 
dominated the ionized volume only at z > 31 and z > 48, respectively. To arrive at 
these results we develop a statistical model for the effect of density correlations and 
discrete sources on reionization and solve it with a Monte Carlo method. 

Key words: galaxies:high-redshift - cosmology:theory - galaxies: formation 



1 INTRODUCTION 

The earliest generations of stars are thought to have trans- 
formed the universe from darkness to light and to have reion- 
ized and heated the intergalactic medium. Knowing how the 
reionization process happened is a primary goal of cosmolo- 
gists, because this would tell us when the early stars formed 
and in what kinds of galaxies. The strong fluctuations in 
the number density of galaxies, driven by large-scale density 
fluctuations in the dark matter, imply that the dense regions 
reionize first, pro ducing on large scales an inside-out reion- 
ization topology l|Barkana fc Loebl |2004 ). This basic pic- 
ture has be en studied and confirm ed with detailed analyti- 
cal models jFurlanetto et alj|2004h . semi- numerical methods 
l|Mesinger fc Furlanettoll2007 ) , and by a variety of large nu- 
merical simulation s (|Mellema et al.ll2006l ; IZahn et aLll2007l ; 
iTrac fc Cer]|2007l ) that solve gravity plus radiative trans- 
fer. The distribution of neutral hydrogen during reioniza- 
tion can in principle be m easured from maps of 21-cm emis- 
sion by neutral hydrogen (jMadau et al.lll997l ). although up- 
coming experiments such as the Murchison Widefield Ar- 
ray (MWAfl and the Low Frequency Array (LOFARfl are 
expected to be able to detect ionization fluctuations only 
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statistically (for revie ws see, e.g., iFurlanetto et at 2006b; 
iBarkana fc Loebll2007l h 

The infancy of cosmic reionization, when only a small 
fraction of the volume of the universe was ionized, is 
of interest for a number of reasons. First, when ioniz- 
ing sources were rare at early times, they are expected 
to have formed separate H II bubbles which if observed 
can be used to study directly th e propertie s of individ- 
ual sources and their surroundings l|Cenll200fJ ). without the 
complications of later times, when overlapping bubbles im- 
ply that galaxy clustering dominates the ionization distri- 
bution and the 21-cm power spectrum. Second, when ion- 
ization fluctuations disappear over much of the universe, 
it becomes possible to use the 21-cm technique for other 
applications including those of fundamental cosmology, 
without the complications of ionization fluctuations which 
are intrinsically non-linear (since the ionization fraction 
varies from to 1). Major such application s include mea- 
surements of the density power spectrum dHogan fc Reesl 
119791 : IScott fc Reeslll99of ), of fluct uations in the Lyct radia- 
tion emitted by the first galaxies ( Barkan a fc Loeb] [ 2005b ; 
Pritchard fc Furlanettol l200d : IChuzhov. Alvarez fc Shapiro] 



20061 ), and of fluctuations in the rate of heating from early 
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X-rays (jPritchard fc Furlanettol l2007l ). If ionization fluctu- 
ations are negligible then the angular anisotropy of the 
21-cm power spectrum makes it possible to measure sepa- 
rately various fluctuation sources, including in particular the 
cosmologically-interesting baryonic density power spectrum 
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( Bark ana fc Loebll2005al ). On small scales, the existence of 
H II bubbles (even when rare) affects the fluctuations in 
Lya and X-ray radiation, producing a small-scale cutoff in 
the 21-cm power spectrum that can be used to detect and 
study the population of ga laxies that formed .just 200 million 
years after the Big Bang (Naoz & Barkana 2008). 

While analytical models and numerical simulations ex- 
ist that can be used to study the later epochs of reioniza- 
tion, the early times are very difficult to investigate. Simu- 
lations, which in general must overcome the huge disparity 
between the large characteristic scales of galaxy clustering 
at high redshift and the s mall scales of individual galax- 
ies (|Barkana fc Loebll2004) . are stretched even further at 
early times, when ionizing sources become very rare and 
even larger cosmological volumes are required in order to 
assemble a reasonable statistical sample. As discussed in de- 
tail below, current anal ytical models based on the model of 
iFurlanetto et all (|2004T l account for galaxy clustering but 
are based on continuous variables and cannot account for 
the fact that galaxies are discrete sources. This discreteness 
becomes a crucial factor in the early stages of reionization, 
when the number of ionizing sources per bubble is small. 
In this limit, Poisson fluctuations also become substantial, 
weakening the correlation between the galaxy distribution 
and the underlying large-scale density fluctuations in the 
dark matter. Discreteness can also play a significant role 
during the central stages of reionization, particularly in the 
case of He reionization by quasars, which are rare sources 
believed to form only in massive halos that correspond to 
many-er density fluctuations at high redshift. These various 
aspects of discrete sources are not accounte d for in current 
analytical models. IFurlanetto fc Ohl ((2008) considered he- 
lium reionization and showed that the continuous models 
break down when discreteness is important. They suggested 
to instead use a pure stochastic Poisson model, without halo 
correlations, when He is less than ~ 50% ionized globally. 

In this paper we develop a model that accounts for dis- 
crete sources as well as density correlations. We solve the 
model with a Monte Carlo method and use it to show that 
galaxy correlations play a major role even in the infancy 
of cosmic reionization. Isolated one-source bubbles do dom- 
inate at sufficiently high redshifts, but the pure stochas- 
tic Poisson model is essentially never a good description of 
the bubble size distribution. In the next section we first re- 
view previous models (section 12. P , then develop ours (sec- 
tion 12. 2p and summarize all the various models whose re- 
sults we later compare (section I2.3J1 . We illustrate our re- 
sults during the infancy of reionization (section 13 - 1 ft and 
then develop an approximate calculation that allows us to 
scan through a wide parameter space of possible reioniza- 
tion scenarios (section |3.2[) . Finally, we illustrate our results 
during later stages of reionization (section 13. 3 p and sum- 
marize our conclusions (section [4} . We assume a standard 
ACDM universe with cosmological parameters that match 
the five-yea r WMAP data and ot her large scale structure ob- 
servations l|Komatsu et al. I l2008h . namely fi m = 0.28 (dark 
matter plus baryons), Qa = 0.72 (cosmological constant), 
Q b = 0.046 (baryons), h = 0.7 (Hubble constant), n = 0.96 
(power spectrum index) and ag — 0.82 (power spectrum 
normalization) . 



2 MODEL 

Analytical approaches to galaxy formation and reionization 
are based on the mathematical problem of random walks 
with barriers. The statistics of a random walk with a barrier 
can be used to calculate various one-point distributions, in- 
cluding the distribution of ioniz ed bubble sizes during reion- 
ization (IFurlanetto et al]|2004h . This distribution indicates 
how likely it is for each scale to determine whether a given 
point is ionized. As such, it indicates the relative importance 
of various scales in reionization, yielding important intuition 
about the internal dynamics of reionization. If bubbles of a 
given radius R are common, this produces a strong corre- 
lation in the neutral fraction (and thus 21-cm emission) on 
a scale ~ R, since the ionization states of two points sep- 
arated by up to R are then often coupled. Calculations of 
the 21-cm correlation function using two-point extensions of 
the mod el yield reasonable agreement with numerical sim- 
ulatio ns l|Furlanetto et alj|2004 IZahn et al.ll2007l ; iBarkanal 
12007ft and indicate that the main feature of the power spec- 
trum during reionization, i.e., enhanced large-scale power, 
indeed appears on scales corresponding to the most likely 
bubble sizes. 

In this section we first review the basic setup of the 
random walk problem in the context of reionization. We 
then show how the standard approach can be generalized 
to solve for the bubble size distribution including Poisson 
fluctuations. 



2.1 Reionization: basic setup 

The basic approa ch for using r andom walks with barriers in 
cosmology follows iBond et all (|l99lf ). wh o used it to rederive 
and e xtend the halo formation model of IPress fc Schechterl 
(1974). In this approach we work with the linear overden- 
sity field 5(x, z) = p{x,z)/p{z) — 1, where x is a comoving 
position in space, z is the cosmological redshift and p is 
the mean value of the mass density p. In the linear regime, 
the overdensity grows in proportion to the linear growth 
factor D(z) (defined relative to z = 0), making it possible 
to extrapolate the initial density field at high redshift to 
the present by multiplication by the relative growth factor. 
Thus, in this paper the density S and related quantities re- 
fer to their values linearly-extrapolated to the present. In 
each application there is in addition a barrier that signifies 
the critical value (as a function of scale) which the linearly- 
extrapolated S must reach in order to achieve some physi- 
cal milestone, which here corresponds to having a sufficient 
number of galaxies within some region in order to fully reion- 
ize it. 

Considering an arbitrary point A in space (at a given 
z), we calculate as follo ws its probability of bei ng inside H II 
bubbles of various sizes (|Furlanetto et al.ll2004r ) . We consider 
the smoothed density around this point, first averaging over 
a large scale or, equivalently, including only small comoving 
wavenumbers k. We then average over smaller scales (i.e., 
include larger k) until we find the largest scale on which 
the averaged overdensity is higher than the barrier; in the 
application to reionization, we then assume that the point 
A belongs to an H II bubble of this size. Mathematically, if 
the initial density field is a Gaussian random field and the 
smoothing is done using sharp fc-space filters, then the value 



of the smoothed 8 undergoes a random walk as the cutoff 
value of k is increased. Instead of using fc, we adopt the 
(linearly-extrapolated) variance S of density fluctuations as 
the independent variable. While the solutions are derived 
in reference to sharp fc-space smoothing, we follow the tra- 
ditional extended Press-Schechter approach and substitute 
real-space quantities in the final formulas. In particular, S 
is calculated as the variance of the mass M enclosed in a 
spatial sphere of comoving radius r. 

The appropriate b arrier for reionization was derived by 
IFurlanetto et al l \2004 ), who noted that the ionized fraction 
in a region is given by x l = £.Fooil , where F co ii is the collapse 
fraction (i.e., the gas fraction in galactic halos) and £ is the 
overall efficiency factor, which is the number of ionizing pho- 
tons that escape from galactic halos per hydrogen atom (or 
ion) contained in these halos, divided by the number of times 
each hydrogen atom in the intergalactic medium must be 
reionized (where this number is assumed to be s patially uni- 
form ) . In the extended Press-Schechter model (|Bond et al.l 
Il99ll ), in a region containing a mass corresponding to vari- 
ance Sr, 

FcoI1 = erfc ( ^ Z) ~ 5 ) , (1) 

\ yj 2(S' m i n — Sr) J 

where Smin is the variance corresponding to the minimum 
mass Mmin of a halo that hosts a galaxy, 8 is the mean den- 
sity fluctuation in the given region, and S c (z) is the crit- 
ical density for halo collapse at z. In reality, the cosmic 
mean halo distribution in simulations is better de scribed by 
the halo mass function o f ISheth fc Tormenl dl999h (with t he 
updated parameters suggested bv lSheth fc Tormenl 1)20021) ). 
However, an exact analytical generalization is not known for 
the biased F co n in regions of various sizes (corresponding to 
Sr) and mean density fluct uations 5. 

iBarkana fc Loebl l|2004h suggested a hybrid prescription 
that adjusts the abundance in vario us regions based o n the 
extended Press-Schechter formula (|Bond et al.l ll991), and 
showed that it fits a broad range of simulation results. In 
general, we denote by f(5 c (z), S) dS the mass fraction con- 
tained at z within halos with mass in the range correspond- 
ing to variance S to S+dS, where S c (z) is the critical density 
for halo collapse at z. Then the biased mass function in a 
region of size R (correspondin g to density variance S r) and 
mean density fluctuation 5 is (|Barkana fc Lo eb 2004) 

hU6c(z),5,S R ,S) = i STi t 5 ; ( t Z ] ,S 2 f PS (S c (z)~8,S-S R ) , 
fps(0c(z),b) 

(2) 

where /ps and /st are, respectively, the Press-Schechter 
and Sheth-Tormen halo mass functions. The value of 
F co n{S c {z), 5, Sr, S) is the integral of / b i as over S, from 
up to the value S m in that corresponds t o the minimu m halo 
mass M m i n or circular velocity V c = y GM mln / R V1I (where 
J?vir is the virial radius of a halo of mass M m i n at z). We 
then numerically find the value of <5 that gives C,F co \\ = 1 at 
each Sr, yielding the exact barrier. Also, in order to com- 
pare with a simpler, analytically-solvable model, we derive 
a linear approximation to the barrier, S(Sr) ~ v + /j.Sr, 
by numerically finding the value of the barrier at Sr — 
and its derivative with respect to Sr. In general, photon 
conservation implies that the mean global ionized fraction 
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should equal x 1 = £/st in terms of the cosmic mean collapse 

frac tion. 

iBarkanal (|2007T ) and IBarkana fc Loebl (|2008h used an 
approximation in which effectively each factor on the right- 
hand side of equation ([2]) was integrated separately over 
S, yielding a simple analytical formula for the effective 
linear barrier. This ap proximation was also assumed by 
IFurlanetto et al. (2006a) when they stated that this hybrid 
prescription does not change the bubble size distribution 
from the pure Press-Schechter case (for a fixed redshift, min- 
imum halo mass, and cosmic mean ionized fraction). Here 
we solve numerically for the barrier using the exact formu- 
las. We show that the previously-used approximation is not 
too accurate, especially at the early stages of reionization 
that are our focus in this paper. 

2.2 The statistics of a random walk with a barrier 
and discrete sources 

The standard approach presented above treats the random 
walks as functions of a continuous variable Sr, and assumes 
a one-to-one correspondence between the value of <5 and the 
ionized fraction x 1 at each scale. The statistical distribution 
of first barrier crossing, which physically corresponds to the 
bubble size distribution, ca n be derived analytical ly for the 
approximate linear barrier l|Furlanetto et aLll2004T ). and for 
the exact barrier can be solved with Monte Carlo simula- 
tions of random wa lks or by solving an integral equation 
jZhang fc Huill2006T l. 

In reality, there are two additional physical constraints 
that are neglected in the standard approach: the ionizing 
sources are discrete, and the ionized fraction (for a given 
value of S in a region) fluctuates due to Poisson fluctua- 
tions in the number of galaxies. The discreteness of ionizing 
sources means that the possible volume of bubbles has a 
minimum value Vbub corresponding to the bubble due to a 
single galaxy hosted by a halo of mass M m i u . Also, the ex- 
pected ionized fraction x % given by the continuous model is 
subject to Poisson fluctuations, as the actual ionized frac- 
tion depends on the number of galaxies. Unlike the standard 
random walk approach, in which the statistics of the walk 
depend only on the barrier expressed as a function S(Sr), 
Poisson fluctuations introduce an explicit dependence on the 
mapping between Sr and scale R. 

In order to include these discrete aspects in the bub- 
ble distribution, we begin with the standard analytical ap- 
proach, which considers the statistics of spherical volumes 
of various sizes R, all about a point A. Given a value 5 on 
a scale R (with corresponding variance Sr), we now treat 
the continuous ionized fraction x l of the previous subsection 
only as an average expectation value. To find its real distri- 
bution, we first calculate the mean expected value (j) of the 
number of ionizing sources within the sphere of radius R. 
This (non-integer) value can be calculated from the integral 
of /bias dS weighted by 1 /M (which yields halos weighted by 
number rather than mass); it depends on z, S m in, 5, and Sr. 
The actual value of j is given by a Poisson distribution with 
mean equal to (j). To find the actual x', we find the mass of 
each of the j halos according to the halo mass distribution 
given by /bias; note that this procedure does not involve a 
single, fixed mass distribution since /bias is a function of S 
and Sr. 
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The complicating factor in this procedure is that we 
cannot treat each scale R independently, since the ionizing 
sources are correlated among the various volumes. This is 
the case first because the densities 8 are correlated, and sec- 
ond because the Poisson fluctuations are correlated, since 
each sphere contains all the galaxies that lie within all 
smaller enclosed spheres. The correlation of the densities 
is dealt with in the standard way reviewed above, where 
small-scale power is added gradually as smaller spheres are 
considered. This makes 5(5*2) dependent on S(Si) if S2 > Si, 
forcing us to start on large scales Sr = and go to smaller 
ones. However, the Poisson fluctuations are correlated in the 
other direction, since a region S2 contains a region Si if 
S 2 < Si. 

The solution is a two-step Monte Carlo method: first, 
we generate the random walk S(Sr), going from Sr = 
to its maximum value (corresponding to the minimum bub- 
ble volume Vbub) in equal steps. At each Sr step, we find 
the mean expected number of galactic halos {j)(Sn) and 
the mean expected total mass of these halos, (Af to t) (Sr). 
Note that the mean expected ionized fraction is {x z )(Sr) = 
C(M to t)(SR)/M(S R ), where M(Sr) is the total mass con- 
tained within the spherical volume of radius R. In the second 
step, we generate the actual ionized fractions starting from 
the smallest scale, Vbub, and working outwards. At Vbub, we 
generate an instance of a Poisson distribution with mean 
(i) (Sr.), yielding an actual integer number j of halos, for 
each of which we find its mass from the appropriate dis- 
tribution of halo number versus mass, derived from /bias- 
Then, for each larger scale (i.e., smaller Sr value), the ad- 
ditional number of galaxies from the last step is on average 
expected to be the difference (Aj) = (j)(Si) — (j)(S2), where 
Si < 5*2 are two consecutive steps in Sr. We find the ac- 
tual difference Aj from a Poisson distribution with a mean 
of (Aj). However, while an actual number of galaxies can- 
not be negative, sometimes the random walk in S gives a 
value (Aj) < 0. In this case we assume that Aj = 0, since 
the number of galaxies already enclosed in a smaller volume 
(corresponding to S2) must also be found in the larger, en- 
closing volume (Si). We do not discard the negative value of 
(Aj) but add it in the next step to the next value of (Aj), 
continuing until we reach a positive expected mean value on 
which we can operate a Poisson distribution. In each step, 
we also keep track of the expected total mass difference, 
(AM) = (Aftot) (Si) - (Mtot)(S 2 ). We slightly modif£| each 
distribution /bias that we use to generate the individual halo 
masses so as to give the correct expected (AM) . This proce- 
dure ensures that on each scale we obtain the correct average 
number of galaxies and correct average galaxy mass, both to 
high accuracy. We note also that in each step in Sr, even if 5 
at the end of the step is below the barrier, there is a chance 
that the random walk hit the barrier during the step. We 
estimate this probability using a linear barrier approxima- 
tion applied separately to each step, and if the walk hit the 
barrier then we raise 5 at the end of the step to the barrier. 
This procedure greatly accelerates the convergence of the 



3 We scale the input M value of the cumulative distribution of 
halo mass M (so that the total probability of having M 
remains unity), with the scaling factor (typically close to unity) 
chosen to yield the correct mean halo mass. 



results as a function of the total number of steps adopted in 
Sr. 

2.3 Summary of models 

We summarize here the various models for the bubble size 
distribution that we consider and compare below. 

(i) Model A: The correct distribution as given by our full 
model. The bubble size distribution is calculated with our 
Monte Carlo method with discrete sources and Poisson fluc- 
tuations, as detailed in section 12.21 We also keep track of 
how many sources are contained in each generated bubble, 
which allows us to find the cumulative volume fraction con- 
tained in bubbles with at least N sources, as a function of 
N. 

(ii) Model B: The exact, continuous barrier (without 
Poisson fluctuations or discreteness). We calculate the non- 
linear barrier S(Sr) numerically, as detailed in section 12.11 
We then derive the bubble size distribution with a Monte 
Carlo method that generates random walks and tracks where 
they first cross the barrier. 

(iii) Model C: A continuous linear barrier approximation. 
We calculate a linear barrier approximation 5(Sr) ~ v+hSr 
numerically, as detailed in section 12.11 W e then derive the 
bubbl e size distribution analytically as in iFurlanetto et all 
l|2004h 

(iv) Model D: The previously-used continuous linear bar- 
rier approximation. Here we apply the additional approxi- 
mation mentioned at the end of section [23] where we noted 
that it gives the same bubble size distribution as in the linear 
barrier approximation of the pure Press-Schechter (rather 
than Sheth-Tormen) model. In this case we calculate ana- 
lytically a linear barrier approximation S(Sr) ~ v + hSr 
and then deri ve the resulting bubble size distribution ana- 
lytically as in IFurlanetto et al.l (12004 ). 

(v) Model E : The pure sto c hastic Poisson model 
suggested by IFurlanetto fc Ohl (120081 ). This model, 
which neglects halo correlations and assumes randomly 
placed, equal-intensity sources, yields an analytical result 
IFurlanetto fc Ohll200ct ) for the cumulative volume fraction 
contained in bubbles with at least N sources. 

Note that the minimum bubble scale is Vbub for models 
A and E, and Vbub/C (corresponding to the scale of the min- 
imum halo mass M m i n ) for models B-D. Also note that we 
have tested our barrier-crossing Monte-Carlo code by com- 
paring it to the analytical solution of a continuous linear 
barrier (Models C and D). We have confirmed precise con- 
vergence, to within a relative error below 1% in the total 
ionization probability, i.e., the total probability of crossing 
the barrier. 



3 RESULTS 

We illustrate our results for a wide range of possible param- 
eters for either hydrogen reionization or helium (full) reion- 
ization. In the latter case, £ is simply interpreted as the over- 
all efficiency factor of producing helium-ionizing photons in 
halos. For hydrogen, minimum halo masses M m i n that are 
often considered are the minimum mass for atomic cooling 
(corresponding to a circular velocity V c ~ 16.5 km/s, where 
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V c — yj GM/R in terms of the halo virial mass and radius) , 
or (at very high redshift) the minimum mass for molecular 
hydrogen cooling (V c ~ 4.5 km/s). Also, much larger min- 
imum masses are possible, V c ~ 35 km/s due to photoion- 
ization feedback (which should affect most of the universe 
by the time reionization is well advanced), or even larger 
values if internal supernova feedback strongly decreases the 
star formation efficiency of low-mass halos at high redshift. 
For helium reionization, assuming it occurs much later, pho- 
toionization feedback affects the source halos from an early 
stage when the density of the assembling matter is still low, 
resulting in a larger V c ~ 80 km/s (i.e., of order the Jeans 
mass). Furthermore, if the observed super-linear local rela- 
tion between halo and black hole mass holds at high redshift, 
then quasars are relatively much brighter in more massive 
halos, increasing the typical halos of helium-ionizing sources 
to V c ~ 200 — 300 km/s. For a given V c , the efficiency ( can 
be chosen to give complete reionizatiorQ x' = f at various 
redshifts z te i. For a fixed z rc i, a larger M m i n implies that 
rarer halos caused reionization, resulting in larger Poisson 
fluctuations. 



3.1 Basic results and comparison with previous 
models 

We begin by considering examples corresponding to an early 
stage (mean ionized fraction x 1 = 1%) of hydrogen or he- 
lium reionization. For hydrogen, we assume atomic cooling 
(V c = 16.5 km/s), with efficiency set to complete reioniza- 
tion (i.e., x % — 1) at z rc i = 7 (implying £ = 19). For helium 
we assume z rc i = 3 and V c = 285 km/s (implying £ = 95), 
which gives the bubbles a minimum size atz~3ofi?=10 
Mpc, about the expected size for the quasars that are ob- 
served to dominate the ionizing photon production at that 
redshift (|Furlanetto fc Oh|[200a) . Figure [T] shows that the 
bubble size distribution obtained from our full model is sub- 
stantially different from the predictions of previous models 
that are based either on a continuous barrier or on a purely 
stochastic Poisson distribution. 

In general these models, based as they are on spherical 
statistics, do not precisely conserve photons, and thus do not 
yield precisely the desired x % = C/st which we have set to 
1%. Indeed, the raw total ionization probability yielded by 
the models is 2.1% (H) and 2.0% (He) for model D, 0.90% 
(H) and 1.4% (He) for model C, and 1.3% (H) and 1.4% 
(He) for model A. Thus, in the figure we compare the rel- 
ative distributions, expressed in terms of the fraction Fi of 
the ionized volume contained in various bubbles. Note that 
model E is defined according to the desired x l , and model B 
(the exact continuous barrier) is mathematically consistent 
in the sense that it yields the correct total x % if the probabil- 
ity is integrated down to V = Vbub/C (We have numerically 
verified this mathematical consistency to a relative error of 
~ 1%). 

Discreteness strongly fails for the continuous barrier 

4 Note that our simple model does not include the likely added 
importance of recombinations towards the end of reionization. 
However, in this paper we do not consider the late stages of reion- 
ization except as a convenient fiducial mark for normalizing £ 
through z re i. 
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V /V bub or N [#sources] 

Figure 1. Cumulative bubble volume distribution as a function 
of V/Vbub! or of the number N of ionizing sources in the bubble. 
Assuming z re ; = 7 and V c = 16.5 km/s for H, and z re ; = 3 and 
V c = 285 km/s for He, we consider x l = 1% (z = 16.3 for H, 6.3 
for He). We compare Fife V), the fraction of the ionized volume 
contained in bubbles with volume ^ V, from our full Monte Carlo 
method with discrete sources (Model A; solid curves) to Fife V) 
from a continuous barrier (Model B; short-dashed curves), a con- 
tinuous linear barrier approximation (Model C; dotted curves), 
and a continuous linear Press-Schechter barrier (Model D; dot- 
dashed curves). We also show Fife N), the fraction of the ionized 
volume contained in bubbles containing TV sources, from our 
full model (Model A; H: squares, He: stars; long-dashed curves 
for N > 10), and from a pure Poisson model (Model E; circles). 

models (both linear and non-linear), in the sense that much 
of the ionized volume in these models is predicted to occur 
inside bubbles below the minimum volume Vbub, especially 
in the case of Helium reionization. Indeed, Fife Vbub) is 
only 55% (H) and 8.2% (He) for model D, 70% (H) and 13% 
(He) for model C, and 82% (H) and 19% (He) for model 
B. Thus, the continuous barrier models fail since they as- 
sign a substantial probability to the unphysical case of frac- 
tional bubbles that are produced by less than one source. Ex- 
pressed differently, the continuous barrier models underpre- 
dict Fife Vbub) since they do not include the Poisson fluc- 
tuations that allow large regions to sometimes reach x % = 1 
even when their mean expected ionized fraction (x % )(Sr) is 
below unity. 

Figure [1] also illustrates the continuous model with a 
linearly approximated barrier, a model use d very commonly 
becau se it yields analytical predictions (|Furlanetto et all 
l2004h . The error of the linear barrier approximation grows at 
small scales, and becomes a 10% error in the barrier height 
at V ~ 0.07Vbub (H) or V ~ 0.02V bu b (He). However, the 
linear barrier approximation becomes relatively accurate on 
scales larger than the scale Vbub corresponding to a one- 
source bubble. On that scale, the height of the linear barrier 
in the examples shown here is only slightly below the height 
of the real barrier (by 2.6% for H and just 0.05% for He), 
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though when i' C 1 the barrier corresponds to a rare ~ 3-cr 
fluctuation on this scale (and rarer still at larger scales), and 
thus small differences in barrier height translate to larger 
differences in Fi. The figure also shows that the pure Press- 
Schechter model (model D) is a rather poor approximation 
to model C. The Sheth-Tormen hybrid model yields more 
large bubbles than the Press-Schechter model, which agrees 
with the expectation based on the Sheth-Tormen mean halo 
mass function, which yields more rare, massive halos than 
does the Press-Schechter mass function. 

While the continuous barrier model extends unphysi- 
cally to V < Vtmb, A does indicate correctly the fact that 
Fife V) declines much more rapidly with V for the He case 
we consider than for H reionization. In fact, we find that 
if we simply cut off the V < Vbub portion and renormalize 
the continuous models relative to V — V bu b (which is not a 
standard way of interpreting these models), then the exact 
and linear barrier models yield nearly identical results, and 
they both yield a reasonable rough estimate to the true bub- 
ble size distribution in the full model. This is illustrated in 
Figure^ which shows the same quantities as in Figure[T]ex- 
cept that all the continuous models have been renormalized 
and are plotted only for V ^ Vbub- For instance, the ratio 
V 1/2 = Fi(V > V bub )/Fi(V > 2Vbub) is 1.18 (H) and 2.33 
(He) in the full model (model A), 1.14 (H) and 2.54 (He) for 
the continuous exact barrier (model B), and 1.15 (H) and 
2.58 (He) for the continuous linear barrier (model C). This 
approach to the continuous models provides a reasonable es- 
timate of the full bubble size distribution in the case of H 
reionization; e.g., Vi/i 00 = Fi(V > V buh )/Fi(V ^ lOOVbub) 
for H is 6.73 in model A, 6.46 in model B, and 6.21 in 
model C, so that here the linear model C, calculated an- 
alytically (i.e., without using Monte Carlo random walks or 
Poisson fluctuations), yields an estimate of Vinoo that is 
within 8% of the true answer according to model A. How- 
ever, for He reionization this approach is much less suc- 
cessful in predicting ratios involving large volumes; e.g., 
V 1/5 = F,(V > V huh )/Fi(V > 5Vbub) for He is 9.9 in model 
A, 16.0 in model B, and 17.6 in model C, and these differ- 
ences increase with V (Figure [2]). 

With our full model (model A), we can also sepa- 
rately predict the distribution by number Fife N). This 
drops more rapidly with N than the distribution by vol- 
ume Fife V) does with V, since large- volume bubbles 
can be produced either by having many sources of mass 
~ M m i n or with a smaller number of individually mas- 
sive halos taken from the high-mass end of the halo mass 
function. Still, Fife N) declines with N much more slowly 
than a pure Poisson model would predict. Indeed, a purely 
stochastic model as suggested by iFurlanetto fc Obi (|200&t ) 
for the early stages of He ionization (or even as late as 
x x ~ 50%), where Poisson fluctuations are assumed that 
are uncorrelated with the underlying density distribution, 
completely fails to descr ibe the results. The ana lytical pre- 
dictions of this model l|Furlanetto fc Oh 2008) yield, for 
x 4 = 1% (for either H or He), Fi(N ^ 2) = 2.0 x 10" 2 
and Fi(N > 3) = 4.4 x 10" 4 (with the latter already out- 
side the range of Figures [1] and [2} . In particular, the ratio 
from the previous paragraph (but applied to the number of 
sources), N 1/2 = Fi(N ^ l)/Fi(N ^ 2), is 1.32 (H) or 3.8 
(He) in the full model, compared to N\/2 = 51 for model E. 
Clearly, density correlations play a substantial role in deter- 
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Figure 2. Cumulative bubble volume distribution as a function 
of V/Vbubi or of the number N of ionizing sources in the bubble. 
Same as Figure [T] except that the continuous models (models B— 
D) have been cut off below V = Vbub and renormalized at that 
point. Note that the curves for models B and C nearly overlap. 



mining the abundance of multi-source bubbles, even early 
on in reionization and even when the process is driven by 
large, rare ionizing sources (such as quasars). 

To help understand the relation of the full model to 
the pure Poisson and to the continuous barrier models, we 
show in Figure the relation between ionization in bubbles 
and the underlying linear density 8. Density fluctuations are 
strongly correlated with ionization, so that the density of 
ionized regions is strongly biased high, and the distribution 
is very different from the standard Gaussian that would be 
expected in a pure Poisson model. However, Poisson fluctu- 
ations allow regions to fully ionize themselves even if their 
density is significantly lower than the barrier, which in a 
continuous model would set the minimum needed 8 for ion- 
ization by internal sources. In particular, the median 8 for 
regions ionized by exactly N sources (where 'exactly' means 
not contained in any larger H II region) represents a fluctua- 
tion of 2.4-ct, 2.57-ct, and 2.61-cr (for H) or 1.9-cr, 2.4-cr, and 
2.7-ct (for He), for N = 1, 2, and 3, respectively. The corre- 
sponding (median) barriers, on the other hand, are 2.907-ct, 
2.909-ct, and 2.922-cr (for H), or 3.2-cr, 3.5-cr, and 3.7-er (for 
He). Thus, the barriers do give a good rough indication in 
each case of whether the 8 distributions for various -/V are 
spaced out or squeezed together. This in turn determines 
whether one-source bubbles are dominant and N > 1 is 
rare, or if multi-source bubbles are at least as common as 
N = 1. 

The continuous model indicates that the main parame- 
ters controlling the relative dominance of single-source bub- 
bles are the effective efficiency £ and the effective slope of 
the power spectrum on the scale i?bub of a one-source bub- 
ble. The efficiency sets the ratio between the scale J?bub of a 
one-source bubble and the scale i? m in from which a galactic 
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Figure 3. Cumulative probability distribution of the linear den- 
sity 5 in units of its standard deviation a = \/S on the relevant 
scale. For either H (dashed curves) or He (solid curves) reioniza- 
tion with the same parameters as in Figure [l] we show 5) for 
regions ionized by exactly one, two, or three sources (from left to 
right in each set of curves). In each case, a circle on the one-bubble 
curve shows the median barrier height on the corresponding scale. 
Also shown for comparison is the cumulative distribution of the 
normal distribution for unconstrained regions (dotted curve). 

halo of mass M m i n was assembled (this ratio equals C^ 3 )- 
Now, the key issue is the relative difficulty of each scale 
achieving self-ionization, when we consider different scales. 
To self-ionize, a region must reach a high enough collapse 
fraction, which according to the extended Press-Schechter 
formula in equation {T}, requires a value of 5 that depends 
on the variance (S m i u — Sr) available for density fluctuations 
inside the regioqj. In order to reach this required value (i.e., 
the barrier), the density has the variance Sr to work with. 
Thus, when we increase the scale (e.g., going from a typical 
one-source bubble to one with two sources), if the fractional 
decline in Sr is more rapid than in (S'min — Sr), then self- 
ionized regions become rarer quickly with increasing scale, 
leading to the dominance of one-source bubbles. This is the 
case when Sr <C S'min, i.e., it requires first that the bubble 
and halo scales differ by a large factor (which requires large 
values of £), and also that the variance depend significantly 
on scale (otherwise, Sr and S m in will be about the same 
even if the corresponding scales are very different). 

More quantitatively, the fractional decline in Sr, over 
the fractional decline in (S m i n — Sr), is (for small changes 
in Sr) equal to (Suiih/Sr — 1). If the power spectrum of 
density fluctuations is approximated as a power law with an 
effective index n over the relevant range of scales, then this 
ratio, which indicates how much harder (in terms of number 

5 The Sheth-Tormen hybrid model alters things slightly, but we 
use the simpler formula here as a rough guide for a qualitative 
understanding. 



of a of the fluctuation) it is to ionize larger scales, is ap- 
proximately £ 1+ ( n ' 3 ) — 1. On small scales, n approaches the 
asymptotic value of —3, making all scales behave roughly 
equally even when ( is relatively large. Note, though, that 
increasing £ increases i?bub and thus brings larger scales into 
play, making the effective n less negative and thus boosting 
the effect of the increased £ on making few-source bubbles 
dominant. This puts a quantitative face on the intuition that 
rare sources tend to create bubbles with small numbers of 
sources. To illustrate, in our H example i? m in = 64 kpc and 
Rbub = 169 kpc, giving n ~ —2.5, while in the He exam- 
ple Rmm = 1.7 Mpc and -Rbub = 7.8 Mpc, giving n ~ —2. 
Thus, He reionization by quasars has both a high efficiency 
and corresponds to a relatively large scale, both of which 
contribute to making small bubbles more dominant, in par- 
ticular the smallest bubbles created by single sources. 

3.2 Approximate calculation 

The current lack of observations at high redshifts leaves ba- 
sic parameters of the galaxy population unconstrained at 
early times. While our model can be used to calculate the 
bubble distribution in any particular case, the need to run 
Monte Carlo trials makes it difficult to explore a large pa- 
rameter space. Thus, an approximate but quick calculation 
is useful for this purpose. In developing such an approxi- 
mation, we focus on determining when one-source bubbles 
dominate the ionizing volume. This can be investigated with 
the ratio N1/2, which is close to unity when multi-source 
bubbles dominate, and is 3> 1 when one-source bubbles do. 
Specifically, this ratio is related to the fraction Fi(N — 1) of 
the ionized volume that is contained in one-source bubbles 
through JVi /2 = 1/[1 - Fi(N = 1)]. 

To construct an approximate calculation of this ratio 
we first adopt the approximation of having equal intensity 
sources, all corresponding to halos having a mass equal to 
the mean expected mass (M). While this approximation 
does not work well for obtaining information on the bubble 
size distribution, we find that it works reasonably for our 
desired ratio involving the distribution of number of sources 
per bubble. We first consider in general the self-ionization 
probability on the scale of a bubble containing j sources 
(with a variance 5* that we approximate as that correspond- 
ing to a volume jVbub), i.e., the probability that a region of 
this size contains at least j sources (regardless of whether or 
not it is contained in some larger bubble). A first attempt 
to calculate this quantity P se if (j) is to calculate the Poisson 
probability of having at least j sources, averaged over the 
normal distribution of 5 on the scale S: 

PmuU) = I d5-^=e- s2 ^ 2S ^P Pois (jx\8, S); > j) , (3) 
J V27tA 

where Pp i s (a; ^ j) denotes the probability of having at 
least j sources in a Poisson distribution with mean a, and 
x % (which also depends on z and S'min) is an approximate 
estimate of (x % ) where we use the same approximation as 
in model D in order to obtain a simple formula. For large 
bubbles, equation [3] for P se if(j) underestimates the self- 
ionization probability, since for a given mean 8 in the re- 
gion, internal density fluctuations increase the variance of 
the number of sources beyond a pure Poisson distribution. 
For j = 2 we can instead calculate a more accurate self- 
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ionization probability by calculating a double integral over 
the joint normal distribution of Si and 62, the mean den- 
sities inside a one-source volume Vbub and inside the sur- 
rounding two-source volume, respectively. Given S± and 82, 
the mean expected number of sources in the two regions 
is ni = x'(z, Smixx, 61, Si) and n 2 = 2x\z, S min , S 2 , S2), re- 
spectively, where Si and S2 are the corresponding variances. 
The probability of self-ionization of the two-source volume 
is then the probability of having at least 2 total sources from 
the sum of a Poisson distribution of mean n\ plus a Poisson 
distribution of mean 712 — ni (except that the latter quantity 
is restricted to be non-negative, a key point which allows the 
larger fluctuations in n\ to contribute). 

Calculating P se if(jO exactly for j > 2 would require at 
least a triple integration, but since j = 1 and j = 2 are most 
important for estimating N 1 / 2 , we simply estimate the self- 
ionization probability for all j > 2 with equation @. Now, 
Pseii(j) for any j is itself only a lower limit for the ioniza- 
tion probability P(N ^ j), since the region may be part of 
a larger H II bubble even if it cannot fully ionize on its own. 
Actually, when one-source bubbles dominate and P(N ^ j) 
drops rapidly with j, regions are much more likely to self- 
ionize than to get ionization help from larger scales, and 
then P S eif(j) becomes an accurate estimate of P(N ^ j). 
However, in order to achieve reasonable accuracy also when 
multi-source bubbles are important, we add a correction to 
each PseifC?) based on the values of P se if(fc) for k > j. In- 
deed, instead of just calculating P se if (fc), which is the prob- 
ability of having at least k sources in a region of size cor- 
responding to k sources, we can separately estimate Pi(k), 
the probability of having exactly I sources in that region, 
using a formula just like equation Q but using the Pois- 
son probability of finding I sources. Then, for any number 
I ^ k sources, we calculate the additional ionization prob- 
ability that was not previously included in P(N ^ j) (for 
each internal volume j < k) using the approximation that 
the I sources are uniformly distributed within the volume 
k. In this way, we estimate the probabilities P(N ^ 1) and 
P(N ^ 2) including the contributions of larger volumes with 
j > 2. When one-source bubbles dominate, higher- j volumes 
have a small effect, but when multi-source bubbles dominate 
the effect adds up, and we cut off j so that P(N ^ 1) does 
not rise above the global ionization fraction x 1 . Actually, 
we find that while the correction from higher- j volumes can 
change each of P(N ^ 1) and P(N ^ 2) by up to a factor of 
a few (giving results much closer to the full model A), the 
relative effect on their ratio is ~ 15% at most. 

Our estimate for N 1/2 is simply P(N ^ 1)/P(N > 2). 
The approximate calculation becomes exact in the limit 
N1/2 — * 00, where our estimated probabilities P S eif(j) be- 
come very small for all j ^ 2, while in the opposite limit, 
when N1/2 — > 1 all quantities become nearly independent 
of j and thus our estimate for the ratio approaches unity, 
also correctly. In practice, from direct comparison with the 
Monte Carlo method at x l ranging from 10 -6 to 1, and at 
ratios A r 1/ / 2 ranging from 1 to 200, we find that our approx- 
imation for this ratio is accurate to ~ 15% (though below 
we extrapolate it beyond the tested range). 

Having developed a quick, relatively accurate calcula- 
tion method, we can use it to explore which areas of pa- 
rameter space will be dominated by one-source bubbles and 
which will form many multi-source bubbles. Figures [4] and 
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Figure 4. Sweep of the parameter space using our approximate 
calculation, showing the relative dominance of one-source com- 
pared to many-source bubbles as indicated by the ratio N 1 ^ 2 = 
P(N ^ l)/P(N ^ 2). For f = 19 or 95, as indicated, we consider 
galactic halos with minimum V c = 4.5, 16.5, 35, 80, or 285 km/s 
(solid curves, from bottom to top). We compare to the case of 
a pure stochastic Poisson distribution (model E; dotted curves). 
Also shown are the locations corresponding to half of the vol- 
ume being in one-source bubbles (horizontal long-dashed line), 
and to 90% in one-source bubbles (horizontal short-dashed line); 
redshifts are indicated at these locations for each case (if it lies 
within the range of the plot). Note also that the various curves 
are not continued below 2 = 3 (for V c = 285 km/s) or z = 6 (for 
the other cases). 

[5]show the ratio Ni/ 2 in the approximate calculation, for x 1 
ranging from 1 down to 10~ 9 , over the whole relevant range 
of source masses, i.e., assuming a minimum V c = 4.5, 16.5, 
35, 80, or 285 km/s, and for four values of the efficiency £, 
19, 95, 580, and 5800. It is interesting to consider the whole 
parameter space, without normalizing to a particular reion- 
ization redshift, since the dominant population of ionizing 
sources at any given redshift may not have similar properties 
to that near the end of reionization, due to the evolution in 
time of chemical, radiative, and hydrodynamical feedbacks. 

The above values of £ are chosen to be particular inter- 
esting, where C, = 19 corresponds to z rc \ = 7 for V c = 16.5 
km/s, and £ = 95 corresponds to z rc i = 3 for V c = 285 km/s, 
which are the H and He reionization examples considered in 
the previous subsection. More generally for star-forming ha- 
los, Population II stars (assumed similar to low-metallicity 
stars forming today) produce ~ 5800 ionizing photons per 
hydrogen atom in stars, while Population III stars (assumed 
to consist of IOOM0, zero- met allicity stars) produce around 
10 times more. Thus, if we assume a maximum star forma- 
tion efficiency of 10% (i.e., that this fraction of the baryons 
in a halo are contained in stars), then if all ionizing pho- 
tons escape out of the dense surroundings of the stars and 
the halo, we get a maximum possible £ = 5800, with Pop III 
stars. A value of £ = 580 can then represent several possibil- 
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Figure 5. Same as Figure[4]but for £ = 580 or 5800, as indicated. 

ities: perhaps only 10% of photons escape, or 100% escape 
but we assume Pop II stars, or we assume Pop III stars 
but with a star formation efficiency of only 1%. The latter 
value is indeed the efficiency expected for the very first, pri- 
mordial Pop III stars in molecular cooling halos. Numerical 
simulations suggest that in each ~ 10 s Mq halo (contain- 
ing, therefore, ~ 1O 4 M0 in baryons), it is likely that only 
a single lOOM© star forms (|Yoshida et al.l [2006h before its 
feedback disrupts the rest of the halo gas and prevents the 
formation of additional stars, at least for some time. Note 
also that these values of £ neglect recombinations, which can 
only lower the effective C, further. 

Figures |4] and [5] imply the general conclusion that ion- 
izing sources produce isolated, single-source bubbles only 
quite early in reionization, when i' < 1. This is a result 
of the fact that while Poisson fluctuations are large when 
we consider just one or two sources, they are strongly mod- 
ulated by halo bias due to the underlying density fluctu- 
ations. Thus, sources are usually found in high-density re- 
gions, which makes it relatively likely to find other sources 
nearby. As sources become rarer at high redshift, the increas- 
ing correlation strength between halos partially compensates 
for the overall low number density of sources, though even- 
tually the sheer rarity of sources does come to dominate. 
As discussed above, increasing V c or ( at a given x r makes 
sources rarer and brings larger scales into play, making it 
easier to form one-source bubbles relative to multi-source 
bubbles. However, only the most extreme case we consider of 
rare, extremely bright sources (V c = 285 km/s and £ = 5800, 
an highly unlikely combination) approaches the results ex- 
pected for a pure stochastic Poisson distribution; the ra- 
tio in the stochastic m odel is A r 1 / 2 = 1/[1 — exp(— 2a; 1 )] 
l|Furlanetto fc Ohll2008T ). 

The Figures also indicate the redshifts when the fraction 
Fi(N = 1) of the ionized volume that is contained in one- 
source bubbles equals 50% (corresponding to N\/2 = 2) or 
90% (N\/2 = 10). In particular, for V c = 16.5 km/s normal- 



ized to produce H reionization at 2 rc i = 7 (i.e., £ = 19), one- 
source bubbles dominate (i.e., Fi(N = 1) > 90%) only above 
z — 57 (outside the plot range), while multi-source bubbles 
become equally important (i.e., Fi(N = 1) = 50%) at red- 
shift 30. Primordial Pop III stars with £ = 580 and V c = 4.5 
km/s also tend to form multi-source bubbles at rather high 
redshifts, with one-source bubbles remaining dominant only 
down to z — 50, and with multi-source bubbles becoming 
equally important at z = 37. On the other hand, for He 
reionization at 2 rc i = 3 with V c = 285 km/s (i.e., £ = 95), 
these milestones are reached at z — 7.7 and z = 5.2, respec- 
tively. Additional cases where these milestones occur outside 
the plot range of the Figures include £ = 95 and V c = 4.5 
km/s, which reaches N1/2 = 10 at z = 62; £ = 19 and 
V c = 35 km/s, which reaches A r 1 / 2 = 10 at z = 37; and the 
faintest example we consider for individual sources, ( — 19 
and V c = 4.5 km/s, which reaches A r 1/ / 2 = 2 at z = 53 and 
does not reach N1/2 ~ 10 eve n at the most lik ely redshift 
(z = 65) of the very first star l|Naoz et alj|2006h . 

If we consider a range of values of £ for halos of a 
given V c , the global ionized fraction x % corresponding to 
a particular milestone (as defined by a particular value of 
Fi(N = 1)) increases with since increasing ( at a fixed 
x l makes sources rarer, while increasing x l (with a fixed £) 
compensates for this by increasing the source number den- 
sity. For each milestone, however, the redshift, which obser- 
vationally is the most directly relevant quantity, behaves in 
a more complicated way, since it is directly related to the 
number density of sources, and thus depends on the ratio 
i l /C- We find that sources with a given V c can only achieve 
a dominance of one-source bubbles at high redshift, almost 
regardless of the efficiency £ (and thus, regardless of the 
reionization redshift). 

Figure [6] shows the minimum z required to achieve var- 
ious values of Fi(N = 1) (as a function of V c ), assuming 
only that the value of £ lies within some wide range. The 
figure shows that while high values of £ do have a larger 
effect on low-V c halos, the minimum redshift is overall rel- 
atively insensitive to the particular range assumed. In par- 
ticular, assuming 10 < £ < 1000, for He reionization by 
quasars (assuming V c ^ 300 km/s), the volume fraction in 
one-source bubbles Fi(N = 1) can be greater than 50% only 
at z > 4.9, 90% at z > 7.3, and 99% at 2 > 9.1. For H reion- 
ization by stars (assuming V c ^ 35 km/s), these milestones 
require 2 > 18, 2 > 23, and 2 > 28, respectively. The gener- 
ation of atomic-cooling halos (V c = 16.5 km/s) can achieve 
Fi(N = 1) > 50% only at 2 > 24, 90% at 2 > 31, and 
99% at 2 > 38. Finally (again assuming 10 < < < 1000), 
the earliest generation of molecular-hydrogen-cooling halos 
(14 = 4.5 km/s) can achieve these milestones only at 2 > 36, 
2 > 48, and 2 > 61, respectively. 

3.3 Later stages 

As reionization advances, eventually the typical bubble size 
encompasses a large number of ionizing sources, reducing 
the importance of discreteness and of Poisson fluctuations. 
Figures[7]and|8]show the cumulative bubble size distribution 
as in Figure [2] but for later stages of reionization. At these 
times, the continuous barrier models still have a significant 
probability at V < Vbub, especially for He reionization by 
quasars; however, if only the V > Vbub portion is considered 
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Figure 6. Minimum redshift (shown in terms of 1 + z) required to 
achieve a dominance of one-source bubbles, for ionizing sources in 
halos with a minimum circular velocity V c (shown over the range 
4.5 — 300 km/s). The minimum redshift shown here is required 
regardless of the reionization redshift or the ionizing efficiency, 
as long as £ is in the range 10—100 (dashed curves), 10-1000 
(solid curves), or 10—10000 (dotted curves). We consider mile- 
stones when the volume fraction in one-source bubbles is 50%, 
90%, or 99% (from bottom to top in each set of curves). 



as in these figures (see also the discussion in section P3 . 1 P . 
then the linear barrier predictions become essentially iden- 
tical to those of the exact barrier, and the predicted bubble 
size distributions of these continuous models are reasonably 
accurate. Specifically, when x 1 = 10%, for the example of 
H reionization, V\i% and Vi/ioo equal 1.08 and 2.28, respec- 
tively, in model A (the full model), 1.06 and 2.18 in model 
B (continuous barrier), and 1.06 and 2.16 in model C (lin- 
ear barrier). For He reionization, VJ/2 and V1/5 are 1.42 and 
2.56 in model A, 1.44 and 2.98 in model B, and 1.44 and 
2.99 in model C. When the universe is 10% ionized, bubbles 
with a small number of sources still play a major role, e.g., 
one and two-source bubbles together account for 19% (H) 
or 60% (He) of the total ionized volume, and the small- iV 
regime is still quite important. 

At the midpoint of global reionization [x 1 — 50%), 
the continuous barrier models approach the full model even 
more (note that the figures at different x % have different y- 
axis ranges). For the example of H reionization, V1/2 and 
^l/ioo equal 1.03 and 1.22, respectively, in model A (the 
full model), 1.01 and 1.19 in model B (continuous barrier), 
and 1.01 and 1.19 in model C (linear barrier). For He reion- 
ization, V1/2 and Vi/ 5 are 1.08 and 1.23 in model A, 1.08 
and 1.24 in model B, and 1.08 and 1.24 in model C. For H 
reionization only 5.2% of the ionized volume lies in one and 
two-source bubbles, but for He this fraction is still 17%. As 
we found in section |3~T1 at x % = 1%, at x 1 = 10% and 50% we 
again see that the pure Press-Schechter model (Model D) is 
a rather poor approximation to model C, and that the pure 
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Figure 7. Cumulative bubble size distribution as a function of 
V/Vb u b, or of the number N of ionizing sources in the bubble. 
Assuming z ro ; = 7 and V c = 16.5 km/s for H, and z re ; = 3 and 
V c = 285 km/s for He, we consider = 10% (z = 12.1 for H, 4.8 
for He). We compare Fife V), the fraction of the ionized volume 
contained in bubbles with volume ^ V, between model A (solid 
curves), model B (short-dashed curves), model C (dotted curves), 
and model D (dot-dashed curves). We also show Fife N), the 
fraction of the ionized volume contained in bubbles containing 
N sources, from model A (H: squares, He: stars; long-dashed 
curves for N > 10) and model E (circles). Note that the curves 
for models B and C essentially overlap. 



Poisson model (Model E) predicts a distribution by number 
Fife N) that falls off much faster with TV than do the true 
distributions (for H or He reionization) according to model 
A. 



4 CONCLUSIONS 

We have developed a model of reionization that adds discrete 
ionizing sources and Poisson flu ctuations to the continuous 
model of lFurlanetto et all (|2004r ). We have shown how to ob- 
tain the distribution of ionized bubbles, versus both bubble 
size and number of ionizing sources, with a two-step Monte 
Carlo method that accounts for both density and Poisson 
correlations among regions of various sizes surrounding a 
given random point in the universe. The bubble size distri- 
bution we obtained differs substantially from previous mod- 
els, but if the continuous barrier model is cut off below Vbub 
(the minimum bubble volume corresponding to a single halo 
of mass M m i n ) then it yields a reasonable rough estimate to 
the true bubble size distribution. More specifically, this esti- 
mate is generally accurate for H reionization even as early as 
a mean ionized fraction x 1 — 1%, while for He reionization 
it works best for small volumes and at later times, and at 
x 1 = 1% is accurate only up to V ~ 3Vbub- Note that with 
the cutoff at Vbub, the linear barrier approximation (which 
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Figure 8. Cumulative bubble size distribution as a function of 
V/Vbubi or of the number N of ionizing sources in the bubble. 
Same as Figure[7] except calculated when x l = 50% (z = 8.7 for 
H, 3.6 for He). Note that the curves for models B and C essentially 
overlap. 



can be calculated analytically) gives nearly identical results 
to the exact continuous barrier. 

Our full model yields a bubble distribution by num- 
ber N that drops more rapidly with N than does the vol- 
ume distribution drop with V , but still, multi-source bub- 
bles are always far more abundant than a pure stochastic 
Poisson model would suggest. This is due to the fact that 
density fluctuations are strongly correlated with ionization 
even when Poisson fluctuations are large. Thus, the density 
of ionized regions is strongly biased high compared to uncon- 
strained regions, but on the other hand, Poisson fluctuations 
allow regions to fully ionize themselves even if their density 
is not as high as would be needed in the continuous barrier 
model. 

The main parameters controlling the relative dominance 
of single-source bubbles are the effective efficiency £ and the 
effective slope n of the power spectrum on the scale of a 
one-source bubble. The ratio of how much harder (in terms 
of number of a of the fluctuation) it is to ionize large bub- 
bles compared to small ones, is approximately proportional 
to £ 1+ ( n / 3 ) — 1. Reionization by rare sources that are mas- 
sive and bright corresponds to having a high £ and to a 
high minimum bubble size, which brings larger scales into 
play, making the effective n less negative and thus making 
it harder to produce multi-source bubbles. 

We have developed a quick, 15% accuracy approximate 
calculation of the ratio N1/2 between the total ionized vol- 
ume and that in multi-source bubbles. This allowed us to 
sweep through the full parameter space of possible halo 
masses and efficiencies of the ionizing sources, and to show 
that sources with a given minimum circular velocity V c can 
only achieve a dominance of one-source bubbles at high 
redshift, regardless of their efficiency or of the reionization 



redshift. In particular, for He reionization by quasars, one- 
source bubbles can dominate (i.e., contain 90% of the ionized 
volume) only at z > 7.3, and fill half the ionized volume at 
z > 4.9, while H reionization by stars can achieve these mile- 
stones only at z > 23 and z > 18, respectively (assuming 
10 < £ < 1000). The generation of atomic-cooling halos can 
place 90% of the ionized volume in isolated bubbles only at 
z > 31 and 50% at z > 24, while the earliest generation of 
molecular-hydrogen-cooling halos can achieve the same only 
at z > 48 and z > 36, respectively. 

We note that reality likely includes even more fluctua- 
tions than included in our Poisson model, since we have still 
assumed that the number of ionizing photons emitted from a 
galactic halo is proportional to its mass. In reality, variations 
in the ionizing efficiency (through spatial or temporal fluc- 
tuations in the star formation efficiency and in the escape 
fraction of ionizing photons), and in the merger histories of 
halos of a given mass (even within a given environment, as 
measured by the average density of a surrounding region) 
will increase the role of (now generalized) Poisson fluctua- 
tions compared to that of galaxy bias due to the underlying 
large-scale density fluctuations. Simple forms of such vari- 
ability can be included in a model of the type that we pre- 
sented, since the ionizing photon outputs from sources are 
added as individual units (which could be generated from 
additional distributions for a given halo mass). In general, 
the model we developed can be used to investigate helium 
reionization and observational prospects for 21-cm observa- 
tions during the infancy of hydrogen reionization. 
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